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Abstract 

We investigate the stability of bounded self-gravitating systems in the canonical en- 
semble by using a thermodynamical approach. Our study extends the earlier work of 
Padmanabhan [ApJ Supp. 71, 651 (1989)] in the microcanonical ensemble. By studying 
the second variations of the free energy, we find that instability sets in precisely at the 
point of minimum temperature in agreement with the theorem of Katz [MNRAS 183, 765 
(1978)]. The perturbation that induces instability at this point is calculated explicitly; 
it has not a "core-halo" structure contrary to what happens in the microcanonical en- 
semble. We also study Jeans type gravitational instability of isothermal gaseous spheres 
described by Navier-Stokes equations. The introduction of a container and the consid- 
eration of an inhomogeneous distribution of matter avoids the Jeans swindle. We show 
analytically the equivalence between dynamical stability and thermodynamical stability 
and the fact that the stability of isothermal gas spheres does not depend on the viscosity. 
This confirms the findings of Semelin et al. [Phys. Rev. D 63 084005 (2001)] who used 
numerical methods or approximations. We also give a simpler derivation of the geometric 
hierarchy of scales inducing instability discovered by these authors. The density profiles 
that trigger these instabilities are calculated analytically; high order modes of instability 
present numerous oscillations that also follow a geometric progression. This suggests that 
the system will fragment in a series of 'clumps' and that these 'clumps' will themselves 
fragment in substructures. The fact that both the domain sizes leading to instability and 
the 'clumps' sizes within a domain follow a geometric progression with the same ratio 
suggests a fractal-like behavior. This gives further support to the interpretation of de 
Vega et al. [Nature, 383, 56 (1996)]. 



1 Introduction 

The thermodynamics of self-gravitating systems is a fascinating subject. It started with 
Antonov (1962) 's discovery that, when a self-gravitating system is confined within a box of ra- 
dius R, no maximum entropy state can exist below a certain critical energy E = — 0.335GM 2 / R. 
This intriguing result was further discussed by Lynden-Bell & Wood (1968) who conjectured 
that for E < — 0.335GM 2 / R the system would collapse and overheat. This is called "gravother- 
mal catastrophe" or "Antonov instability" . Lynden-Bell & Wood have related this phenomenon 
to the very particular property of self-gravitating systems to possess negative specific heats. The 
gravothermal catastrophe picture has been confirmed by sophisticated numerical simulations 
(Larson 1970, Cohn 1980, Lynden-Bell & Eggleton 1980) and is expected to play a crucial role 
in the evolution of globular clusters. It is found that the collapse proceeds self-similarly (with 
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power law behaviors) and that the central density becomes infinite in a finite time. This insta- 
bility has been known as "core collapse" and many globular clusters have probably experienced 
core collapse (Binney & Tremaine 1987). In the case of dense clusters of compact stars (neutron 
stars or stellar mass black holes), the gravothermal catastrophe can lead to the formation of 
supermassive black holes of the right size to explain quasars and AGNs (Shapiro & Teukol- 
sky 1995). Statistical mechanics is also relevant for collisionless stellar systems (e.g., elliptical 
galaxies, dark matter,...) undergoing a violent relaxation (Lynden-Bell 1967, Chavanis et al. 
1996, Chavanis & Sommeria 1998, Chavanis 1998a,2001a). In particular, the inner regions of 
elliptical galaxies are close to isothermal and this is an important ingredient to understand de 
Vaucouleurs' R 1/A law (Hjorth & Madsen 1993). 

On a theoretical point of view, the stability of isothermal spheres has been first investigated 
by Katz (1978) with a very powerful method extending Poincare's theory of linear series of 
equilibrium. He found that instability sets in precisely at the point of minimum energy. This 
stability analysis was reconsidered by Padmanabhan (1989) who studied the sign of the second 
variations of entropy and reduced the problem of stability to an eigenvalue equation. This leads 
to the same stability limit as Katz but the method of Padmanabhan provides in addition the 
form of the perturbation that induces instability at the critical point. It is found that this 
perturbation presents a "core-halo" structure. 

The analysis of Padmanabhan (1989) was performed in the microcanonical ensemble in 
which the energy is fixed. The microcanonical ensemble is probably the most relevant for 
studying stellar systems like elliptical galaxies or globular clusters (Binney & Tremaine 1987). 
Indeed, apart from a slow evaporation, these systems can be assumed isolated so the evolution 
conserves energy E and mass M . In addition, from the viewpoint of statistical mechanics, 
only the microcanonical ensemble is rigorously justified for non extensive systems, as discussed 
in the review of Padmanabhan (1990). However, it is always possible to define formally a 
canonical ensemble or a grand canonical ensemble in which the temperature is fixed instead of 
the energy. As suggested by de Vega et al. (1996a, 1996b) these ensembles may be suitable 
for describing the cold interstellar medium where the temperature is imposed by the cosmic 
background radiation at T ~ 3K in the outer parts of galaxies, devoid of any star and heating 
sources (Pfenninger et al. 1994, Pfenninger & Combes 1994). In particular, by working out 
the statistical mechanics of the self-gravitating gas, de Vega et al. (1996a, 1996b) have shown 
that self-gravity can provide a dynamical mechanism to produce the fractal structure of the 
interstellar medium. They used the same approach to explain the fractal structure of the 
universe (de Vega et al. 1998), assuming that galaxies have reached a quasi-thermodynamical 
equilibrium like in the early work of Saslaw & Hamilton (1984). 

For self-gravitating systems, it is well known that the thermodynamical ensembles do not 
coincide in the whole range of parameters (Padmanabhan 1990). Using toy models, Lynden- 
Bell & Lynden-Bell (1977) and Padmanabhan (1990) demonstrated that the region of negative 
specific heats allowed in the microcanonical ensemble is replaced by a phase transition in the 
canonical ensemble. This phase transition separates a dilute "gaseous" phase from a dense 
"collapsed" phase. Since these toy models are not very realistic, the self-gravitating gas was 
also studied in a meanfield approximation. In this viewpoint, an isothermal sphere is stable if 
and only if it is a local maximum of an appropriate thermodynamical potential (the entropy in 
the microcanonical ensemble and the free energy in the canonical ensemble) . As expected from 
physical grounds, phase transition occurs when the gaseous sphere ceases to be a local maxi- 
mum of this potential and becomes a saddle point. In this paper, we investigate the stability of 
isothermal gaseous spheres in the canonical ensemble by studying the sign of the second varia- 
tions of the free energy. Our analysis is a direct extension of Padmanabhan (1989)'s approach 
in the microcanonical ensemble. The two studies therefore provide a unified description of the 
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stability of isothermal spheres in the meanfield approximation in terms of thermodynamical 
potentials. 

For a long time, the thermodynamics of self-gravitating systems was considered exclusively 
in a meanfield approach (or with toy models). However, recently, de Vega & Sanchez (2001a,b) 
have developed a rigorous statistical mechanics of self-gravitating systems by using field theo- 
retical methods. In particular, they evidenced the existence of a thermodynamic limit in which 
the number of particles iV and the volume R 3 go to infinity keeping N/R fixed. This very 
unusual thermodynamic limit proves to be appropriate to non extensive systems. By using 
Monte Carlo simulations and analytical calculations, they showed that the meanfield approx- 
imation correctly describes the thermodynamic limit except near the critical points where a 
phase transition occurs. They also derived the equation of state of a self-gravitating gas instead 
of assuming it, as is done in the meanfield treatments. Therefore, their work fully justifies the 
studies of previous authors and specifies their range of validity. 

This paper is organized as follows. In section ^ we introduce a meanfield description of 
the system in the canonical ensemble and show that critical points of free energy J at fixed 
temperature T and mass M correspond to isothermal spheres like those studied in the context 
of stellar structure (Chandraskhar 1942). We show that there is no global maximum of free 
energy. There is not even a local maximum of free energy in an unbounded domain: unbounded 
isothermal spheres have an infinite mass! We restrict therefore our analysis to the case of self- 
gravitating systems confined within a spherical box of radius R. In this case, there exists 
local maxima of J (metastable states) if the normalized temperature 77 = ^ is less than 
2.52 and the density contrast 1Z = p(0)/p(R) < 32.1. Critical points of free energy with 
density contrast TZ > 32.1 are unstable saddle points. For 77 > 2.52, there are not even critical 
points of free energy: in that case, the system is expected to undergo a phase transition and 
collapse. This "isothermal collapse" is the counterpart of the "gravothermal catastrophe" in 
the microcanonical ensemble (see Figs. [l|-@). 

In section ^, we study the sign of the second variations of free energy by using the methods 
of Padmanabhan (1989) introduced in the microcanonical ensemble. We show analytically that 
instability sets in precisely at the point of minimum temperature in agreement with the theorem 
of Katz (1978). The perturbation that induces instability at this point is calculated explicitly; 
it has not a "core-halo" structure contrary to what happens in the microcanonical ensemble. 

In section |, we study Jeans type gravitational instability of isothermal gaseous spheres 
described by Navier-Stokes equations. The introduction of a container removes the problems 
associated with an infinite homogeneous medium and avoids the Jeans swindle. We show 
analytically the equivalence between dynamical stability and thermodynamical stability and the 
fact that the stability of isothermal gas spheres does not depend on the viscosity. This confirms 
the findings of Semelin et al. (2001) who used numerical methods or approximations. We also 
give a simpler derivation of the geometric hierarchy of scales inducing instability discovered by 
these authors using sophisticated renormalization group technics (Semelin et al. 1999). This 
provides a more illuminating interpretation of their results. 

In section [5], we make speculations about the fragmentation and the fractal structure of an 
isothermal self-gravitating gas. We distinguish between the Jeans length Lj defined with the 
mean density and the King's length Lx defined with the central density. If we fix the Jeans 

length, instability occurs for R > ^/ 2 ' 5 3 2 '" Lj and is marked by the absence of critical point 

of free energy (i.e., hydrostatic equilibrium) above this threshhold. In that case the system 
is expected to collapse without fragmenting. If we fix the King's length (or core size), a first 
instability occurs for R = 8,9 3 9,,, Lx- Above this threshold critical points of free energy still 
exist but they are unstable saddle points. Secondary instabilities occur at larger box radii 
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that asymptotically follow a geometric progression R n ~ [10.74...] n L^. The density profiles 
that trigger these high order modes of instability are calculated analytically. They present 
more and more oscillations whose zeros also follow a geometric progression r n ~ [10. 74. ..]"/?. 
The profile that destabilizes the singular isothermal sphere has an infinite number of nodes! 
Each oscillation can be interpreted as a 'germ' in the langage of phase transition and the above 
picture suggests that the system will fragment into a series of 'clumps'. It is expected that these 
'clumps' will evolve by achieving higher and higher density contrasts, and finally fragment in 
turn into substructures. This yields a hierarchy of structures fitting one into each other in a 
self-similar way. This picture is given further support by the fact that both the domain sizes 
inducing instability and the zeros of the perturbation profile in each domain follow a geometric 
progression with the same ratio. This double-geometric progression may explain in a natural 
way the fractal structure of a self-gravitating gas like the interstellar medium and the large 
scale structures of the universe. This gives further support to the interpretation of de Vega et 
al. (1996a, 1996b, 1998) who first realized the importance played by self-gravity in building a 
fractal distribution of matter. 



2 Thermo dynamical stability of self-gravitating systems 
in the canonical ensemble 

2.1 The free energy 

Consider a system of N particles, each of mass m, interacting via Newtonian gravity. The 
particles can be galaxies, stars, atoms, etc... We assume that the system is non rotating and 
non expanding. Let /(r, v, t) denote the distribution function of the system, i.e. /(r, v, t)d 3 rd 3 v 
gives the mass of particles whose position and velocity are in the cell (r, v; r + d 3 r, v + d 3 v) at 
time t. The integral of / over the velocity determines the spatial density 



f<t\. 



On the other hand, in the meanfield approximation, the total mass and the total energy of the 
system can be expressed as 

(2) M = Nm = / pd 3 r, 



(3) E = \j ^ 2rf3rrf3v + \ J = K + W , 

where K is the kinetic energy and W the potential energy. The gravitational potential <£> is 
related to the star density by the Newton-Poisson equation 

(4) A$ = 4nGp. 

Finally, the Boltzmann entropy is given by the standard formula 

(5) S = -k [ — In ^-d 3 rd 3 v, 

J m m 

which can be obtained by counting the number of microstates corresponding to a given macrostate 
and taking the logarithm of this number (Ogorodnikov 1965). 
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We shall work in the canonical ensemble in which the temperature is fixed, allowing the en- 
ergy to fluctuate. In that case, the relevant thermodynamical potential is the Massieu function 
related to the Helmholtz free energy F = E — TS by J = —^F. Hence 

(6) J = S-±E. 

At equilibrium, the system is expected to be in a state that maximizes the Massieu function 
d^) for a fixed total mass M (in the following we shall call J the free energy, although this is 
only the free energy up to a negative proportionality factor). 

2.2 The isothermal gaseous spheres 

Following Padmanabhan (1989) 's procedure, we start to maximize the free energy J for a given 
density field p(r). This yields the Maxwell-Boltzmann distribution 

(7) f={^w) 

which is a global maximum of J with the previous constraint. Substituting this optimal dis- 
tribution function in Eqs. ©-(i), we can express the energy and the entropy in terms of the 
spatial density in the form 

(8) E = ^NkT + ± J p$d 3 r, 



(9) 



S 3N 3N , 

- = 1 In 

k 2 2 



2M 



m 



^ln^i 3 r. 

m m 



We can now determine the variations of J around a given density profile p(r). To second order 
in the expansion we get 



m 



5J=-— I 6p[ ln-^ + 1 )cfr- — 



m 



(10) 



T 



m 
1 



2p 



- / ^5pd 6 r - — / 5p5$d a T. 



2T 



Introducing a Lagrange multiplier a to satisfy the conservation of mass, the condition that 
J is an extremum is written (to first order) 



= 5 J - a5M 



d 6 r 



$ kin 
- + — [ In — 

m 



a 



Sp. 



T m 

This condition must be satisfied for any variations Sp. This yields the Boltzmann distribution 
(12) p = Ae-P*, 

where we have set 



(13) 



P 



m 

~kT' 



Therefore, the Boltzmann distribution fll2|) is a critical point of free energy. This does not 
insure, however, that it is a maximum of J. It is not even clear that the extremum problem 
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leading to Eq. ([TJ) has a solution. Indeed, the gravitational potential that appears in Eq. (|T2| ) 
must be determined self-consistently by solving the mean field equation 

(14) A$ = A-nGAe-^, 

obtained by substituting the density fli~2| ) in the Poisson equation (|j), and relating the Lagrange 
multiplier A to the constraint M. As we shall see, this problem does not always have a solution. 
When a solution exists, we must consider the sign of the second order variations 5 2 J to determine 
whether it is a maximum or not. 



2.3 Absence of global maximum of free energy 

It is easy to show that there is no global maximum of free energy. To prove this result, we 
just need to consider a homogeneous sphere of mass M and radius R. The total energy of this 
sphere is 

and its free energy 

(16) J = -Nk\n[ -iVATn M . 

K ; 2 \ m J J 5TR 

From the above formula it is clear that J goes to +oo when we spread the density to infinity 
(R — > +oo) or when we contract the system to a point (R — > 0). Since the mass is conserved 
in this process no global maximum of free energy can exist. 

In addition, there is not even a local maximum of free energy in an unbounded domain: 
unbounded isothermal spheres have an infinite mass. For non rotating systems, the equilibrium 
states are expected to be spherically symmetric. In that case, the Boltzmann-Poisson equation 
flT4]) can be rewritten 



(17) LM r ^) =i7rGAe -^. 

r z dr \ dr J 

This equation has been studied extensively in the context of isothermal gaseous spheres (Chan- 
drasekhar 1942). We can check by direct substitution that the distribution 

(18) $ s (r) = hn(27rG(3Ar 2 ), p s (r) ' 



(3 y r " 2-nGfir 2 ' 

is an exact solution of Eq. (|T7| ) known as the singular isothermal sphere (Binney & Tremaine 
1987). Since p ~ r -2 at large distances, the total mass of the system M = J + °° p Anr 2 dr is 
infinite! More generally, we can show that any solution of the meanfield equation ( |lTD behaves 
like the singular sphere as r — > +oo and has therefore an infinite mass (Chandrasekhar 1942). 

We shall avoid the infinite mass problem by confining artificially the system within a spheri- 
cal box of radius R. It is only under this simplifying assumption that a rigorous thermodynamics 
of self-gravitating systems can be carried out. Physically, the maximum cut-off determines the 
scale at which other processes intervene to limitate the spatial extent of the system. It is clear 
from Eq. (|l^) that the introduction of an upper cut-off does not remove the absence of global 
maximum of free energy (the system can always increase J by collapsing to a point). However, 
the presence of a wall at radius R allows the existence of local maxima of free energy (metastable 
states). 
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2.4 The equilibrium phase diagram 

The micro canonical and canonical ensembles yield the same critical points, i.e. the critical 
points of entropy at fixed mass and energy and the critical points of free energy at fixed mass 
and temperature coincide. Only the onset of instability will differ from one ensemble to the 
other. The thermodynamical parameters for bounded isothermal spheres in the meanfield 
approximation have been calculated by Lynden-Bell & Wood (1968) and we shall directly use 
their results. 
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Figure 1: Stability diagram for isothermal spheres in the microcanonical ensemble (fixed E). 



To that purpose, we introduce the function ip = /3($ — $ ) where <3> is the gravitational 
potential at r = 0. Then, the density field can be written 



(19) 



P = Poe 



where po is the central density. Introducing the notation £ = (4:7iGPpoY^ 2 r, the Boltzmann- 
Poisson equation (|T7| ) reduces to the standard Emden form (Chandrasekhar 1942) 



(20) 



Eq. (|20[ ) has a simple analytic solution, the singular sphere 

(21) e"* = 1. 

The regular solutions of Eq. (p0|) must satisfy the boundary conditions 

(22) V = $ = 0, 

when £ = 0. These solutions must be computed numerically. However, their asymptotic 
behaviors are well-known 



(23) 
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Figure 2: Stability diagram for isothermal spheres in the canonical ensemble (fixed T). 



(24) 



2 r a 



cos 



-oo 



The curve intersects the singular solution (^T|) infinitely often at points that asymptotically 
increase geometrically in the ratio 1 : e 27T ^ = 1 : 10.74.... This property will have important 
physical consequences in section |5|. 

In the case of bounded isothermal systems, the solutions of Eq. (|20|) are terminated by the 
box at different radii given by a = (4irGpp ) 1/2 R. Lynden-Bell & Wood (1968) show that the 
reduced temperature and reduced energy can be expressed in terms of a by 



(25) 



0GM 
R 



at/)' (a), 



(26) 



A 



ER 



GM 2 2at/)'(a) i//(a)< 



For each value of the normalized inverse temperature ^ , we can solve Eq. (|25|) to get 
a. Substituting the result in Eq. (|26|), we deduce the corresponding value of the normalized 
energy —-M3i- We can then determine the equilibrium diagram E — (3 (Figs. [l|-|2]). The critical 



points of entropy or free energy form a spiral. This curve is parametrized by a that goes from 
to +oo as we spiral inward. Instead of the parameter a, it may be more relevant to introduce 
the density contrast 



(27) 



n 



Po 
p(R) 



which gives a more physical parametrization of the solutions. For a — > 0, we can use the 
asymptotic behavior ip ~ |a 2 of the potential at the origin [Eq. (^)] and we find ~ ^ — > 
0, —-Mk ~ — 2^2 ~^ —°° an d TZ —> 1. This corresponds to high temperatures T — > +oo. In this 



case, self-gravity is negligible with respect to thermal motion and the system behaves like an 
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ordinary gas with uniform density. We check, by eliminating a from the two foregoing relations, 
that E = %M-. This is indeed the equation of state for an ideal gas without self-gravity. When 
we proceed along the spiral, the density contrast increases monotonically and goes to +00 at 
the center of the spiral. Using the asymptotic expansion e~^^ ~ \ of the potential at large 
distances [Eq. (||)], we find ^ -»■ 2, -> \ and K ~ ^ -> +00. We can check by a 

direct calculation that the ending point (1/4, 2) of the spiral corresponds to the singular sphere 

m 



2.5 The Milne variables 

It will be convenient in the following to introduce the Milne variables (u, v) defined by (Chan- 
drasekhar 1942): 

(28) u = ^— and v = 

ip' 

Taking the logarithmic derivative of u and v with respect to £ and using Eq. fl2"0|), we get 

1 du 1 . 

(29) __ = . (8 -„-„), 

(30) i* -^-l,. 

Taking the ratio of the foregoing equations, we find that the variables u and v are related to 
each other by a first order differential equation 

. . udv u — 1 

(31) 



v du u + v — 3 

Therefore, by using the Milne variables, the degree of the meanfield equation (|20|) has been 
reduced from two to one. As discussed extensively by Chandrasekhar (1942), this property 
is related to the homology invariance of the solutions of the Emden equation. The solution 
curve in the (u, v ) plane is well known and is ploted in Fig. |^. The curve is parametrized by 
£. Starting from the point (u, v) = (3, 0) corresponding to £ = 0, the solution curve spirals 
indefinitely around the point (u s ,v s ) = (1,2), corresponding to the singular sphere, as £ tends 
to infinity. All isothermal spheres must necessarily lie on this curve. For bounded isothermal 
spheres, £ must be terminated at the box radius a. Clearly, the spiral behavior of the (w, v) 
curve (and also the curve of Figs. |l|-^D can be ascribed to the oscillating behavior of the solution 
($Z§ as f -> +00. 

It turns out that the normalized temperature and energy can be expressed very simply in 
terms of the values of u and v at the normalized box radius a. Indeed, writing uq = u(a) and 



Vq = v(a) and using Eqs. (E^)-(B6|), we get 



(32) A = -(l~u 

v V 2 



(33) rj = v . 



The relation ([32|) was previously noticed by Padmanabhan (1989). The intersection between 
the lines defined by Eqs. fl32|)-(|33"D and the spiral in the (u, v) plane determines the value of a 
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Figure 3: The (u,v) plane. All isothermal spheres must necessarily lie on the spiral. There 
exists solutions only for r\ < 2.52 and A < 0.335. 

corresponding to a given energy or temperature. As noticed by Padmanabhan, for A > A c = 
0.335 there is no intersection. Thus, in the microcanonical ensemble, no isothermal sphere can 
exist if —-§^2 > 0.335. This result can also be read from Fig. |l]and was discovered by Antonov 
(1962). For sufficiently small energies, the system is expected to collapse and overheat in 
agreement with the "gravothermal catastrophe" picture (Lynden-Bell & Wood 1968). Similarly, 
considering Eq. (^), we find that there is no intersection for > v max = 2.52 (where v max 
corresponds to u — 1). This result can also be read from Fig. |2|. In the canonical ensemble, a 
gaseous sphere is expected to collapse below a critical temperature kT c = c £^r (Lynden-Bell 
& Wood 1968). 

2.6 The stability analysis of Katz 

More precisely, it is possible to prove the following stability result: when a self-gravitating 
system is confined within a box and maintained at a constant temperature T, local maxima 
of free energy exist only for ^ M < 2.52; they have a density contrast 1Z = p(0)/p(R) < 32.2 
(and a < 8.99). Critical points of free energy with density contrast 1Z > 32.2 (and a > 8.99) 
are unstable saddle points. For > 2.52, there are no critical points of free energy. 

The stability of the solutions can be deduced from the topology of the E — [3 curve by using 
the method of Katz (1978) who has extended Poincare's theory of linear series of equilibrium. 
The parameter conjugate to the free energy with respect to the inverse temperature (3 is — E = 
(jm)M,R- Then, if we plot E as a function of (3 (we just need to rotate Fig. || by 90°) we have 
the following results: (i) a change of stability can occur only at a limit point where (3 is an 
extremum (dE/d(3 infinite) (ii) a mode of stability is lost when the curve rotates clockwise and 
gained otherwise. Now, we know that for T sufficiently large the solutions are stable because, 
in this limit, self-gravity is negligible and the system behaves like an ordinary gas. From 
point (i) we conclude that the solutions with 1Z = p(0)/p(R) < 32.2 are stable. As the curve 
spirals inward for TZ > 32.2, more and more modes of stability are lost. In this respect, the 
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singular sphere, at the end of the spiral, is the most unstable solution []. It can be noted that 
instability sets in precisely when the specific heat C = dE/dT becomes negative (in the range 
32.2 < 1Z < 709). General considerations indicate that negative specific heats are forbidden in 
the canonical ensemble (Padmanabhan 1990). By contrast, thermally isolated systems can have 
negative specific heats: in the microcanonical ensemble, the solutions with 32.2 < 1Z < 709 
are stable (Katz 1978, Padmanabhan 1989). This clearly demonstrates that the statistical 
ensembles do not coincide for self-gravitating systems: the region of negative specific heats in 
the microcanonical ensemble is replaced by a phase transition (an "isothermal collapse" ) in the 
canonical ensemble (Padmanabhan 1990). It is at the verge of this phase transition that the 
meanfield approximation ceases to be valid as demonstrated by de Vega et al. (2001a, 2001b). 



2.7 Analogy with critical phenomena 



In this section, we study some analogies with the theory of critical phenomena. Further analo- 
gies are discussed in Chavanis et al. (2001) where a simple dynamical model for studying phase 
transitions in self-gravitating systems is proposed. 
The specific heat can be written 



(34) 

Using Eqs. 
(35) 



c _dE_ _ _k_p2<^_ _ N^ 2 — 
dT m d(3 dr\ 



)(]33|) (ID©, we easily get 

dq_ _ vo 
da a 



(uo ~ 1), 



(36) 



dA 

da 



2avo 



(Aul + 2u v - 11m + 3) 



Therefore, the specific heat is determined as a function of a by the expression 



(37) 



C = Nk 



4-Uq + 2uqVq — llwo + 3 
2{u - 1) ' 



Of course, C —>■ oo if drj = and C —>■ if dA = 0. On the other hand, for the ideal gas 
(uo,v ) = (3,0), we recover the well-known result C = ^Nk. The specific heat is ploted as a 
function of a on Fig. ^. 

At the critical point rj = r] c at which C diverges, we write a = a c + e and expand r\ and A 
in terms of e, using Eqs. (|55|)-(|3"6j) and u Q (a c ) = 1, v (a c ) = r\ c . To lowest order, we have 



(38) 



V - r] c 



d 2 i] 
da 2 



Vc 

2a 2 c 



(2-?? c 



. . (dA\ ri c -2 

39 A _ A = e = ± e, 

\daj ac a c r] c 

where A — 0.197 is the value of A at a = a c . Eliminating e from these expressions, we find 
the relation between r] and A close to the critical point 

(40) lfc _„ = _!l_ ( A_A 0)2 . 



1 Another application of Katz theorem where the spiral unwinds and the stability is regained is given by 
Chavanis & Sommeria (1998) for self-gravitating fermions. 
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Figure 4: Evolution of the specific heat as a function of a. The specific heat diverges for the 
first time at a c = 8.99 and is negative in the range 8.99 < a < 34.4. For a = 0, C = \Nk 
(perfect gas). 




Figure 5: Evolution of the specific heat as a function of A near the point A . The specific heat 
diverges like C ~ —(A — Ao)^ 1 . It is negative for A > A and tends to \Nk for A — > — oo 
(high energies). These results are supported by the numerical simulations of Cerruti-Sola et al. 
(2001) in the microcanonical ensemble who solved the gravitational A^-body problem in a finite 
container. 
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Figure 6: Evolution of the specific heat as a function of rj near the critical point r\ c . The specific 
heat diverges like C ~ ±(^ c — i])^ 1 ^ 2 . It tends to \Nk for 77 — > (high temperatures). The 
region of negative specific heat (lower curve) is forbidden in the canonical ensemble. 



Then, using Eq. (|34]), we find explicitely 



:v^(a-a )-\ 

1/2 



(41) 



c 



±Nk 



ri c (ri c -2) 
2 



(rj c -rj) 1/2 . 



The specific heat diverges as a function of the energy with an exponent v = 1 and as a 
function of the temperature with an exponent v' = 1/2 (see Figs. |||-|6]). It becomes negative for 
Ao < A < A c . This region of negative specific heats is allowed in the micro canonical ensemble 
but not in the canonical ensemble. 

It is also of interest to determine the behavior of the central density po near the critical 
point (C. Sire, private communication). The central density is related to a by the relation 
a = (AttGPpoY^R or > alternatively, 

(42) Po 



To first order in e, we have 
(43) 



Po ~ Po 



2a c M 



AnR 3 r] c 

where Pq = a 2 M/ (47ri? 3 r/ c ) is the value of the central density at the critical point r] c . Using 
Eqs. (H)-H), we get 

(44) p ° p0 = ) r / Mi/2 



Pi 



± 



1 - i 

Vc 



Quite remarkably, this stationary result can also be derived from dynamical considerations (see 
Chavanis et al. 2001). Note also that, using Eqs. (^) (|28|) , we find 1Z C = ot 2 /rj c1 which implies 
p c (R) = M/(4vr J R 3 ). 
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3 An extension of Padmanabhan's stability analysis 
3.1 The condition of thermodymamical stability 

The stability of the gaseous spheres can be studied by extending Padmanabhan (1989) 's stability 
analysis to the canonical ensemble. This will provide the profiles of the perturbations that 
induce instability. 

A critical point of free energy is a local maximum if, and only if, the second variations 

,2 T _ k f (6p) 2 3 1 f « 



(45) ^ J= __y^_ A __y ff ^, 

are negative for any variation 5p that conserves mass to first order. This conservation law 
imposes 

(46) J 5pd 3 r = 0. 

We shall restrict ourselves to spherically symmetric perturbations since only spherically sym- 
metric perturbations can induce instability for non rotating systems. Following Padmanabhan 
(1989), we introduce the variable q defined by 

^ ^ ^ Anr 2 dr 

Physically, q represents the mass perturbation q(r) = 5M{r) = J Q r Anr' 2 5p(r')dr' within the 
sphere of radius r. The boundary conditions on q are thus 

(48) q(0) = q(R) = 0. 
Substituting the foregoing expression for q(r) in Eq. (f45|), we obtain 

(49) ^=-^ri**-Af * (fw, 

2T J dr m J 8npr 2 \dr J 



Integrating by parts and using the boundary conditions (|4q) , we get 

-2 t 1 f R J 5 ^^ k f R d ( 1 dq 



^ ^ 2T J Q ^ dr m J ^ dr \87r pr 2 dr ) 

Then, using Gauss theorem in the form 

cW$ Gq 



2 



(51) 
we find 



(52) 6J =^ ~h dr + - Qt 5 2Tl dr > 

21 J r z m J dr \6ir pr z dr ' 



dr r 



,2 j 1 f R Gq 2 j_ , k f R d ( 1 dq 



or, equivalently, 



(53) 5 J = - / drq 

2 Jo 



1 



G k d ( 1 d 
+ 



Tr 2 m dr \ATipr 2 dr 



14 



The second variations of free energy can be positive (implying instability) only if the differential 
operator which occurs in the integral has positive eigenvalues. We need therefore to consider 
the eigenvalue problem 



(54) 



k d / 1 d \ G 



mdr\ Anpr 2 dr J Tr 



q x (r) = Xq x (r), 



with <7a(0) = q\(R) = 0. If all the eigenvalues A are negative, then the critical point is a 
maximum of free energy. If at least one eigenvalue is positive, the critical point is an unstable 
saddle point. The point of marginal stability rj = rj c is determined by the condition that the 
largest eigenvalue is equal to zero. We thus have to solve the differential equation 

(55) k d ( 1 dF ) I GF ^ - 

m dr \47r pr 2 dr J Tr 2 

with F(0) = F(R) = 0. Note that this equation corresponds to that found by Padmanabhan 
(1989) in the microcanonical ensemble with V = 0. 

3.2 The point of marginal stability 

Introducing the dimensionless variables of section |2.4j , we can rewrite the foregoing equation in 
the form 

(56) ±(e_^dF\ m_ Q 

with F(0) = F(a) = 0. As shown by Padmanabhan (1989), this equation can be solved without 
solving explicitly the Emden equation (|20|) . Let us introduce the differential operator 



r d fe^ d\ 1 

(57) c ^{eTi) + e- 



Using Eq. (p0|) , it is readily established that 

(58) «€V) - + * - ^* x e-) + 4 = «. 
Similarly, we find 

(59) = ^-(3 - «.') + (e~* = - f + fe"* = if. 



These identities suggest to seek a solution of Eq. flBq) in the form 

(60) F{£) = Cl e 3 e^ + c 2 e^' 

where c\ and c 2 are constants. Substituting the foregoing expression for F(£) in the differential 
equation (|56|) , we find that c\ + c 2 = 0. The function F(£) can therefore be reexpressed as 

(61) F(e) = d(£ 8 e-*-eV)- 
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The boundary condition F(0) = is automatically satisfied (we can show furthermore that Eq. 
( |5"T| ) is the only solution consistent with this boundary condition). The condition F(a) = will 
determine the point on the spiral of Fig. ^ at which the gaseous spheres become unstable (i.e., 
the point of marginal stability). From Eq. (63), we get 



(62) 



a 3 e-^ - 



a 2 ijj'(a) = 0. 



This can be expressed in terms of the Milne variables fl28|) as 



(63) 



M = 1. 



Considering Eq. (|31|) , we see that the condition (|63| ) determines the points for which v is 
extremum (dv/du = 0). Since r\ = t> according to Eq. (|33|), we deduce that a change of stability 
occurs each time that the temperature is extremum (see Fig. |7]) in complete agreement with 
the theorem of Katz (1978). In particular, the first instability (for which the largest eigenvalue 
A becomes positive) corresponds to r\ c = v max = 2.52 (for this value a = 8.99). We show in 
the next section that the secondary instabilities occur at values of a that follow asymptotically 
a geometric progression. The profiles of the perturbations that induce these instabilities are 



discussed in section 3.4 




Figure 7: Evolution of the inverse temperature i] as a function of the scaled radius a. A mode 
of stability is lost each time that rj is extremum. This figure is the counterpart of Fig. 3 of 
Padmanabhan (1989) in the microcanonical ensemble. 



3.3 Secondary instabilities 

A mode of stability is lost each time that uq = u(a) = 1. Now, it is easy to have an analytical 
estimate of m(£) for large values of £. Introducing the asymptotic behavior ( p4|) of ip in the 
Milne variable u defined by Eq. ( p8|) we get for £ — > +oo: 



(64) u 



l + ^cos(^lne + 5) 



1 + [v / 7sin(^ lnf + 5) + cos(^ lnf + 5)} ' 
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Therefore, the condition uq — 1 corresponds to values of a satisfying 



fy/7, \ 3 



(65) tanl-^- In a + 5j = — j= (a — > +oo), 
or equivalent ly 

(66) — In a n + 5 = arctan(3 / v7) + ^tt (n integer) . 

This determines a succession of values that follow the geometric progression 

(67) a n ~ = [10.74...]" (n -> +oo, integer). 

At these points a new eigenvalue becomes positive implying secondary instabilities. Clearly, this 
geometric progression can be traced back to the 'curious' asymptotic behavior of the solutions 
governing the isothermal gas sphere (see Eq. (124]) ) which intersects the singular solution (|2l|) 
infinitely often at points that asymptotically increase geometrically in the ratio 1 : 10.74. 

A mode of stability is lost each time that a achieves one of the values given by Eq. fl6~7P . If 
one prefers, the same criterion can be expressed in terms of the density contrast TZ. Using Eq. 
(|27| ) and the asymptotic expansion fl24|) we get 1Z = ~ ^- for a — > +oo. Hence, the result 
is translated in 



(68) Tl n ~eV7 = [115.5.. .f (n -> +oo, integer). 

3.4 Profiles of the perturbation at the critical points 

In this section, we study the form of the perturbations that trigger instability at the critical 
points. Following the method of Padmanabhan (1989), it is possible to describe the behavior 
of these perturbation profiles without numerical work. Indeed, according to Eq. (H7p, the 
eigenfunction associated with the eigenvalue A = can be written 

< 69 > - = TH§ 

Po 4vr£ 2 d£ 

where F(£) is given by Eq. (|5TD . Substituting this result in Eq. ( |B9"D and simplifying the 
derivative with the aid of Eq. (|20|) we obtain 



(70) ^ = £(2.^). 



Introducing the Milne variable v defined by Eq. (p8|) , we can rewrite the foregoing equation in 
the form 

The density perturbation 5p becomes zero at the point(s) & such that v(£i) = 2. The 
number of zeros is therefore given by the number of intersections between the spiral in the 
(u,v) plane and the line v = 2 (see Fig. ||). For the first mode of instability («i = 8.99), 
we need to follow the spiral until the first extremum of v, which corresponds to the point at 
which the box terminates (recall that u(a\) = 1, see section |3]^). In that case there is only 
one intersection £i with the line v = 2. Therefore, the first mode of instability does not present 
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Figure 10: Second mode of instability corresponding to «2 = 98.0. 



a "core-halo" structure (Fig. |9|), contrary to what happens in the microcanonical ensemble 
(Padmanabhan 1989). For the second mode of instability (a^ = 98.0), we need to follow the 
spiral until the second extremum of v. The intersection with the line v = 2 now determines 
two zeros £i and £2, the first one being the same as for the first mode of instability. The profile 
of the second mode of instability therefore has a "core-halo" structure (Fig. [10]). When we 
spiral inward, more and more intersections are obtained. Therefore the high order modes of 
instability present more and more oscillations. For these modes, it is easy to determine the 
asymptotic positions of their zeros. From Eqs. ( p8|) and fl24]) , we have for £ — ► +00: 

( 72 ) v = 2+ 2^ X 



\fl sin {^~7^~ In £ + + cos 



2 



Thus, substituting v = 2 in Eq. ([72]), we get 

(73) — ln£j + 5 = — arctan(l/V / 7) + ire (i integer). 
Therefore, the zeros follow asymptotically the geometric progression Q 

(74) & ~ = [10.74..]* (i -> +00, i integer). 

Note that the perturbation that destabilizes the singular sphere (at the end of the spiral) has 
an infinite number of oscillations! 



4 Dynamical stability of isothermal gaseous spheres 

We shall now investigate the dynamical stability of isothermal spheres described by Navier- 
Stokes equations and compare the results with the thermodynamical approach. This problem 

2 It might be recalled, parenthetically, that the positions of the planets in the solar system also follow a 
geometric progression, but with a different ratio q ~ 2 (Titius-Bode law). This is recognized to be an effect of 
scale invariance (Graner & Dubrulle 1994, Dubrulle & Graner 1994). There is, however, no real connexion with 
the present work and planets may have formed quite differently (see Chavanis (2000) and references therein). 
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was previously addressed by Yabushita (1968) in the isobaric ensemble and by Semelin et al. 
(2001) in the canonical ensemble. We shall prove analytically and with no approximation the 
equivalence between dynamical and thermodynamical stability and the fact that the stability 
of isothermal spheres does not depend on the viscosity. We shall also describe the behavior of 
the velocity profiles at the critical points. 

The equations of the problem are the equation of continuity, the equation of motion and 
the Poisson equation 

(75) ^ + V(pv) = 0, 

(76) — + (vV)v = —Vp - V$ + ^Av + -(C + ?)V(Vv), 
at p p p 3 

(77) A$ = AnGp. 

They must be completed by an equation of state that we take of the form p = p—T where T is 
constant. We also recall that we work within a box of radius R. For astrophysical fluids, the 
Reynolds numbers are so huge that we can neglect the viscosity in the Navier-Stokes equation. 
We shall take therefore r\ = ( = in most of the calculations but we will also prove that 
viscosity does not change the onset of instability. 

Clearly, the stationary solutions of Eqs. (|75|) (|76|) (f77|) correspond to the isothermal gaseous 
spheres studied in the previous sections. It should be recalled that the condition of mechanical 
equilibrium corresponds to the balance between the pressure force and the gravitational force. 
For an isothermal gas, this yields 

/ % ^dp rf$ 

We now consider a small perturbation around a stationary solution and write 

(79) v = 6v(r,t), p = p + 5p(r,t), 

(80) p = p + 5p(r,t), $ = $ + <5$(r,t), 

where the bar refers to the stationary solution (in the following we shall drop the bar). The 
linearized equations for the perturbations are 

rMv kT 

(81) p-— = V5p - <5pV$ - pV5$, 

at m 

(82) ^ + V(p5v) = 0, 

(83) A5$ = AirG5p. 

Since we work with exact stationary solutions of the fluid equations, that are inhomogeneous, 
we do not have to advocate the Jeans swindle (Binney & Tremaine 1987). In this sense, our 
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approach is more rigorous than the classical treatment of Jeans starting from an infinite and 
homogeneous medium Q. 

In the following, we restrict ourselves to radial perturbations. It is known that non-radial 
perturbations do not lead to new instabilities. Writing the time dependance of the perturbation 
in the form Sv ~ e xt , 5p ~ e xt ,..., the equations of the problem become 



(84) 



, r kTdbp r d$ d5<$> 
XpSv = 5 p — - p-—. 

m dr dr dr 



(85) 



\Sp+ \^-{pr 2 5v) = 0, 



(86) 



1 d 

r 2 dr 



,d5<$> 
dr 



AnG5p. 



Introducing the notation ([47D and using the Gauss theorem (^l|) we see that the Poisson equation 



^6| ) is automatically satisfied. The continuity equation flSB]) becomes 

A 



»7) 



dq ^ 1 d , 2 ^ , 
47rr 2 dr r 2 dr 



0. 



This equation is readily integrated. Using the boundary condition q{0) 

A 



0, we get 



(88) 



5v 



Airpr 2 



Substituting this result back into Eq. (84]), we obtain 



(89) 



A 2 



kT d 



47rr 2 ^ m dr 



1 dq 



+ 



dq d$ Gp 



Anr 2 dr J Airr 2 dr dr 

Using the condition of mechanical equilibrium (|78|) , the foregoing equation can be rewritten 

k d ( 1 dq\ Gq A 2 



(90) 



mdr \4npr 2 dr J Tr 2 AnpTr 



This is similar to the eigenvalue equation (Bl) associated with the second variations of the free 
energy. In particular they coincide for marginal stability A = 0. Therefore, dynamical and 
thermodynamical instability occur at the same point in the series of equilibrium. This was not 
obvious a priori since the Navier-Stokes equations without viscosity conserve the free energy. 
Taking into account a finite viscosity, we obtain instead of Eq. , 

k d / 1 dq\ ^ Gq A 2 ^ (C + ^) ^ 

m dr \4:iipr 2 dr J Tr 2 AnpTr 2 pT 3 dr \_r 2 dr ^4^ 



1 d ( Xq 



For A = the viscous term cancels out. Therefore, viscosity does not change the onset of 
instability, nor does it alter the profile of the perturbation that triggers the instability. 

The profiles of the density perturbation have already been discussed in section |3.4j . The 
profiles of the velocity perturbation at the critical points are given by Eq. fl88f ). In order to 
work with dimensionless variables, we introduce the speed of sound 

dp k 



(92) 



dp 



-T 



rn 



P' 



3 It should be recalled, however, that the Jeans procedure is justified in a cosmological context if we take into 
account the expansion of the universe. 



21 




5 

Figure 11: Velocity profile for the first mode of instability corresponding to a.\ = 8.99. 




Figure 12: Velocity profile for the second mode of instability corresponding to a 2 = 98.0. We 
clearly see the separation between a core (Sv < 0) and a halo (5v > 0). 
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Note that c s can be written c s = \L K /t dyn where td yn = (47rGp ) 1//2 is a time scale comparable 
with the dynamical time (Binney & Tremaine 1987). Then 

(93) - = -^>F(0, 
where we have set A' = \td yn - Using Eqs. (j61"|)-([28"|), we find 

(94) — = -—c^'e^iu- 1). 

C S 4:71 

We assume that we are just at the onset of the instability (A' = + ) so that Eq. ([94]) is 
applicable with A' > (the velocity profile at marginal stability, i.e. A' = 0, is simply 8v — 0). 
Since ip'(0) = 0, the velocity always vanishes at the center of the sphere 5v (£ = 0) = 0. The 
other zero(s) are determined by the condition u(£j) = 1. Since the critical points of instability 
are also characterized by u(a) = 1, see section |37|, we deduce that the velocity perturbation 
always vanishes at the wall: 5v (£ = a) = 0. This result is to be expected since the mass does 
not leave the sphere. The other zeros are determined by the intersections between the spiral in 
the (u, v) plane and the line u = 1 (see Fig. |8|). They correspond precisely to the values of a at 
which a new mode of instability occurs, i.e. £j = c^. It is therefore straightforward to determine 
the number of zeros in the velocity profile. For the first instability there are only two zeros at 
£ = and £ = ot\. For the second instability there are three zeros at £ = 0, £ = a± and £ = 
For the n th instability there are n + 1 zeros £ = 0, £ = cki, £ = a 2 ,..., £ = a n . The exact velocity 
profiles given by Eq. fl94|) are displayed on Figs. [11-12 for n = 1, 2. They are in agreement with 



the numerical results of Semelin et al. (2001). For high order modes of instability, the zeros 
follow asymptotically the geometric progression (p?!). Returning to dimensional variables, the 
zeros of the n th mode of instability are located at r = and = (aii/at n )R (i = 1, ...,n). 



5 Fragmentation and fractal structure of an isothermal 
self-gravitating gas 

We shall now express the previous results in a more physical form and make some speculations 
about the fragmentation and the fractal structure of an isothermal self-gravitating gas. 

5.1 The King's length 

Introducing the King's radius (Binney & Tremaine 1987) 



(95) 



J K 



9 ^ 1/2 



the parameter a defined in section |27| can be rewritten 

(96) a = 3^-. 

L K 

The King's radius gives a good estimate of the core radius of an isothermal sphere. In the 
following, we shall consider that the core radius Lk is fixed and that the domain size R is 
increased. Clearly, this amounts to increasing the parameter a along the spiral. From the 
study of section EO, we find that the onset of instability corresponds to 



8 99 

(97) R > ——Lk- 
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For larger radii there are still critical points of free energy but they are more and more unstable. 
According to Eq. flSTD, a new mode of instability occurs at radii that follow asymptotically the 
geometric progression 

(98) R n ~ [10.74... } n L K . 

This result is clearly related to that of Semelin et al. (1999) who first noticed the appearance 
of a hierarchy of scales connected to the instability of isothermal gaseous spheres. However, 
the connection with their work is not straightforward. Semelin et al. (1999) approximate the 
regular solutions of Eq. ( |TTD by the analytical profile ( |i~8D and introduce an arbitrary short 
cut-off distance 5 to avoid the central singularity. By contrast, there is no cut-off at small 
scales in our theory since we consider regular spheres that are exact solutions of the Emden 
equation (although we do not have analytical expressions for them). In this sense, there are no 
approximation nor indetermination in our theory. In fact, an "effective" cut-off is played by 
the King's length (or core radius) but it is not arbitrary and depends on physical parameters, 
the central density po and the temperature T. For r ^> Lk, the density profile behaves like 
r~ 2 so that regular gaseous spheres can indeed be approximated by a r~ 2 profile plus a cut-off 
5 at small scales as done by Semelin et al. (1999). Our approach suggests to take 5 as the 
King's length. We have therefore given a new derivation of their result which is considerably 
simpler and does not introduce sophisticated renormalization group technics. It also avoids the 
introduction of an arbitrary cut-off. This provides, hopefully, an easier interpretation of their 
intriguing results. 

According to section if we increase the domain size along the series of equilibrium (i.e. 
maintaining a fixed core radius) smaller and smaller regions become unstable. The modes 
that trigger these instabilities present numerous oscillations whose zeros also follow a geometric 
progression. These oscillations can be regarded as sort of 'germs' leading to the formation of 
'clumps' preceding the fragmentation of the isothermal gas. It is expected that these 'clumps' 
will evolve by achieving higher and higher density contrasts, and finally fragment in turn into 
substructures. This yields a hierarchy of structures fitting one into each other in a self-similar 
way. This picture is given further support by the fact that both the domain sizes that induce 
instability and the zeros of the perturbation in a given domain follow a geometric progression 
with the same ratio (see Eqs. fl9"8"D([71p). This property may explain in a quite natural fashion 
the fractal structure that a self-gravitating medium (e.g., the interstellar medium, the large 
scale structures of the universe,...) can build under certain conditions. 

The hierarchy of scales that appears in an isothermal gas is due intrinsically to the oscillating 
nature of the equilibrium phase diagram (see Figs. |l|@,0). This curious behavior has been known 
for a long time in the context of stellar structure (Chandrasekhar 1942) but it had not been 
apparently related to the geometric progression of scales inducing instability. 



5.2 The Jeans length 

We can also introduce another length scale defined with the average density p = J^j- of the 
system instead of the central density. This is the Jeans length (Binney & Tremaine 1987) 

/ 9 \ 1/2 

(99) L ^te) ' 

The reduced temperature r\ can be expressed in terms of the Jeans length as (de Vega & Sanchez 
2001a) 

, . (3GM (R 

(100) , = <1_ = 3(- 
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In the following, we shall consider that the Jeans length Lj is fixed and that the domain size R 
is increased. Clearly, this amounts to increasing the parameter t] along the vertical axis of Fig. 
H Using the results of section p76| , we find that the gaseous sphere becomes unstable when 



(101) 



R > 



2.52.. 



For larger radii the system will collapse because there are no critical point of free energy (i.e. 
hydrostatic equilibrium). In that case there are clearly no secondary instabilities. Indeed, when 
the condition ( J10l| ) is realized we leave the spiral and loose the geometric hierarchy of scales 
associated with its turning points. Since the mode of instability at the critical point does not 
oscillate (see section jOl) , we expect that the system will collapse without fragmenting, probably 
like in the study of Penston (1969). 



5.3 Fragmentation or not fragmentation? 

The fragmentation of the isothermal self-gravitating gas and the fractal structure that it can 
generate were previously discussed by de Vega et al. (1996a, 1996b, 1998) and we here improve 
and complete their physical arguments. In particular, the analytical calculation of the modes 
of instablity (without approximation), the geometric progression of both the domain sizes and 
the 'clumps' sizes, and the precise physical picture that emerges from these results are new. In 
addition, we give arguments explaining why the instability of a self-gravitating gas does not 
always lead to a fractal structure (indeed globular clusters or elliptical galaxies, for example, do 
not show a hierarchical structure). In our model, the selection of the regime is related to how 
the size of the domain compares with the Jeans length or the King's length. If R > 0.916...Lj, 
we expect a collapse without fragmentation but if R > 2.996. ..Lk (with a sufficiently large 
ratio) the medium is expected to fragment into several 'clumps'. 
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Figure 13: Evolution of the dimensionless number M as a function of the scaled radius a. This 
number quantifies the inhomogeneity of an isothermal self-gravitating gas. 

The distinction between the Jeans scale and the King's scale is therefore essential in the 
present context and had not been noticed previously. In particular it is not the Jeans lengths 
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but the King's lengths that follow a geometric progression [see formula ([98D 1 contrary to the 
claim of Semelin et al. (1999). These two scales only coincide for a perfect gas with no self- 
gravity that is homogeneous. Indeed, for a — > 0, we have Lj ~ Lk ~ — — > +00. By contrast, 

for a — > +00 (corresponding to the singular sphere) we get L K ~ ^ — > and Lj — > \J\R (we 

have used the asymptotic expansions of section |2.4|) . To appreciate the difference between the 
two length scales, it is convenient to introduce a dimensionless number 

(102) M=^. 

From the above formulae, we have the equivalent expressions 

1/2 / \ 1/2 

Po \ a I a ■ 



\ P J y/OT] \6ip'{a 

For a — > (perfect gas), A/" ~ 1 and for a — > +00 (singular sphere), A/" ~ ^ — > +00. Therefore 
the number A/" quantifies the inhomogeneity of the self-gravitating gas. The evolution of N(a) 
along the spiral is ploted on Fig. [H| From Eq. ( |103| ), Af is just the ratio between the central 
density and the average density to the power 1/2. 



6 Conclusion 

In the first part of the paper, we have investigated the stability of isothermal spheres in the 
canonical ensemble using a thermodynamical approach. This is a natural extension of the 
work of Padmanabhan (1989) in the microcanonical ensemble. As is well known, there is 
no hydrostatic equilibrium for an isothermal gas below a critical temperature kT c = 
(Lynden-Bell & Wood 1968). This corresponds to a situation in which the domain size becomes 
comparable with the Jeans length, i.e. R > 0.916...Lj. In that case, we expect the system to 
collapse. By analyzing the second variations of free energy, we have found that the point 
of minimum temperature is precisely the point at which the free energy ceases to be a local 
maximum and becomes a saddle point in the series of equilibrium. This vindicates the theorem 
of Katz (1978) and provides in addition the form of the perturbation that induces instability. 
Contrary to the microcanonical ensemble, it has not a "core-halo" structure, nor multiple 
oscillations. We expect therefore the system to collapse without fragmenting. On the other 
hand, if we start from isothermal gaseous spheres with high density contrasts 1Z > 32.2 (this 
corresponds to R > 2.996. ..L K ), such spheres are unstable saddle points and probably fragment 
into small 'clumps' (although they would be said to be stable by a naive Jeans argument since 
R < 0.916...Lj). The number of clumps, presumably associated with the number of zeros 
in the mode of instability, increases with the density contrast or, equivalently, as the domain 
size becomes larger and larger with respect to the King's radius. The distinction between the 
Jeans length and the King's length is essential in the present context and had not been noticed 
previously. The Jeans length is related to the parameter r\ determining whether an hydrostatic 
equilibrium state exists or not (for an isothermal sphere) while the King's length is related to 
the parameter a parametrizing the series of equilibrium. They are very different parameters. 

In the second part of the paper, we have extended Jeans instability criterion to the case of 
inhomogeneous bounded gaseous spheres described by Navier-Stokes equations. This treatment 
avoids the well-known difficulties associated with an infinite and homogeneous medium (Jeans 
swindle). Quite remarkably, we have provided an elegant analytical solution of this more general 
problem without making any approximation. 
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Now, the question that naturally emerges is what happens to a gravitationally unstable 
gaseous configuration. The collapse of isothermal gas spheres has been investigated by Penston 
(1969) using Navier-Stokes equations (without viscosity). He found an analytical solution 
for which the collapse is self-similar and develops a finite time singularity (i.e., the central 
density becomes infinite in a finite time). This only partially answers the above question. 
Indeed, his solution is just a particular solution of the fluid equations and other solutions (e.g., 
fragmentation) may be possible. In particular, it is not clear whether his solution is stable with 
respect to non spherically symmetric perturbations. In addition, Penston (1969) works in an 
infinite domain for which the isothermal spheres have infinite mass and are always unstable. 
Therefore, he could not evidence a phase transition between an equilibrium "gaseous" state 
and a "collapsed" state, depending on the value of the temperature. 

The thermodynamics of self-gravitating systems is so intriguing that we have been tempted 
to explore this problem anew from a dynamical point of view (Chavanis et al. 2001). We 
have introduced a numerical algorithm in the form of a relaxation equation (Chavanis 1996, 
2001a; Chavanis et al. 1996) constructed so as to increase monotonically entropy at fixed 
energy in the microcanonical ensemble and free energy at fixed temperature in the canonical 
ensemble. With this numerical algorithm, we have covered the whole bifurcation diagram in 
parameter space and we have checked, by an independant method, the stability limits of Katz 
(1978). We have verified that the density profile that induces instability at the critical point 
has a "core-halo" structure in the microcanonical ensemble but not in the canonical ensemble, 
in agreement with the study of Padmanabhan (1989) and the present study. We have also 
checked that the number of oscillations in the perturbation profile increases when we start 
from unstable isothermal spheres with high density contrasts. With this algorithm, we have 
explored the "bassin of attraction" of the stable isothermal spheres. Since they only are local 
maxima of entropy or free energy, an initial condition can either relax towards a metastable 
equilibrium state or collapse, depending on its topology. Hence, the control parameter A or r\ is 
not sufficient to characterize completely the final state of the system even in the case when an 
equilibrium exist. This depends whether the initial condition lies in the "bassin of attraction" 
of the equilibrium solution or not. The complete characterization of this bassin of attraction is 
difficult but we have given some examples. With our numerical algorithm, it will be possible 
to include rotation in order to compute entropy maxima that are not spherically symmetric 
and determine the influence of angular momentum conservation on the collapse of a rotating 
self-gravitating system (Chavanis & Rosier, in preparation). Some interesting works have been 
done in that direction (Lagoute & Longaretti 1996, Laliena 1999, Fliegans & Gross 2001 and 
Lynden-Bell 2000) but the complete description of rotating isothermal gaseous spheres is still 
in its infancy stage. 

It is also of first interest to test numerically the speculations that we have made concerning 
the development of the instability, i.e. if the system will collapse as a whole or if it will fragment 
and break into substructures. These two regimes are clearly separated, at least in the theoretical 
model that we have considered here (Antonov model). We have suggested that if we start from 
a solution of the Emden equation with a high density contrast, it should break in a series of 
'clumps' associated with the oscillations of the density perturbation that we have calculated in 
the linear regime. Obviously, if we want to follow the nonlinear evolution correctly, we must 
relax the hypothesis of spherical symmetry. 

In this paper, we have assumed that the particles can be treated by classical mechanics and 
that they are point-like. The extension of our work to general relativity (for, e.g., neutron stars) 
is almost straightforward (Chavanis 2001b). We find that there is no hydrostatic equilibrium if 
the system size is smaller than a multiple of the Schwarzschild radius. We also show analytically 
that the point of smallest radius corresponds precisely to the point at which the series of 
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equilibrium becomes unstable. Secondary modes of instability are found; they follow a geometric 
progression like in the Newtonian case but with a different ratio. When the speed of light goes 
to +00, we recover the results of the present study. It is seen that relativistic effects favour 
gravitational instability, i.e. instability occurs sooner than in the Newtonian limit. On the 
other hand, the case of self-gravitating fermions enclosed within a box has been considered 
by Chavanis & Sommeria (1998). Their study was formulated in the context of the "violent 
relaxation" of collisionless stellar systems introduced by Lynden-Bell (1967) but the results also 
apply to quantum particles such as massive neutrinos in Dark Matter models (Ingrosso et al. 
1992) or degenerate electrons in White Dwarf stars (Chandrasekhar 1959). When degeneracy 
is accounted for, there exists a global entropy maximum for all values of energy (Chavanis & 
Sommeria 1998, Robert 1998). It has a "core-halo" structure with a degenerate core and a 
dilute atmosphere. Depending on the degeneracy parameter, there can be a "gravothermal 
catastrophe" at E = -0.335GM 2 /R but the core ceases to shrink when it becomes degenerate. 
Considering a classical gas with a short distance cut-off a essentially leads to the same results 
(see, e.g., Padmanabhan 1990, Follana & Laliena 2000). 

For astrophysical purposes, it is still a matter of debate to decide whether collisionless 
stellar systems like elliptical galaxies are degenerate (in the sense of Lynden-Bell) or not. Since 
degeneracy can stabilize the system without changing its overall structure at large radii, we 
have suggested that degeneracy could play a role in galactic nuclei (Chavanis & Sommeria 
1998). The recent simulations of Leeuwin and Athanassoula (2000) and the theoretical model 
of Stiavelli (1998) seem to go in that direction: if the nucleus of elliptical galaxies contains 
a (primordial) black hole, degeneracy must be taken into account and can explain the cusps 
observed in the center of galaxies. This form of degeneracy is also relevant for massive neutrinos 
in Dark Matter models where it competes with quantum degeneracy (Kull et al. 1996). 

All together, these results suggest that the statistical mechanics and the thermodynamics 
of self-gravitating systems can have quite relevant astrophysical applications at different scales 
of the universe: stars (Chandraskhar 1959), interstellar medium (de Vega et al. 1996a, 1996b), 
globular clusters (Lynden-Bell & Wood 1968), elliptical galaxies (Lynden-Bell 1967, Hjorth 

6 Madsen 1993, Chavanis & Sommeria 1998), dark matter (Ingrosso et al. 1992, Kull et al. 
1996), cosmology (Saslaw & Hamilton 1984, de Vega et al. 1998). ..The same ideas of statistical 
mechanics have been introduced in two-dimensional turbulence to explain the formation and 
maintenance of coherent vortices, like Jupiter's Great Red Spot, which are common features of 
large-scale geophysical and astrophysical flows (see Bouchet & Sommeria 2000, Bouchet et al. 
2001 and references therein). The analogy between the statistical mechanics of two-dimensional 
vortices and stellar systems has been described extensively by Chavanis (1996, 1998b, 2001c). 
This analogy concerns not only the equilibrium state (the formation of large-scale structures) 
but also the relaxation towards equilibrium (Chavanis et al. 1996, Chavanis 2001d) and the 
statistics of fluctuations (Chavanis & Sire 2000). These ideas may also have applications in 
the context of planet formation where large-scale vortices present in the solar nebula could 
efficiently trap dust particles to form the planetesimals and the planets (Barge & Sommeria 
1995, Tanga et al. 1996, Bracco et al. 1999, Chavanis 2000). Accordingly, the statistical 
mechanics of these nonlinear media seems to be able to account for the fascinating process of 
self-organization in nature. 
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